library(lmtest)
library(sandwich)
library(foreign)
library(survey)



# Load data ---------------------------------------------------------------

svyz <- read.csv("../data/asylum_data.csv")


# Weight pre-processing: Make all weights 1 -------------------------------

svyww <- subset(svyz,!is.na(svyz$weight))
svyww$weight <- 1



# Country vectors ---------------------------------------------------------

countries<-c("Germany","Hungary","Sweden","Austria","Norway",
             "Switzerland","Denmark","Netherlands","Greece","Czech Republic","Italy",
             "Poland","Spain","France","United Kingdom")

countries2<-c("Germany","Hungary","Sweden","Austria","Norway",
              "Switzerland","Denmark","Netherlands","Greece","Czech Republic","Italy",
              "Poland","Spain","France","United Kingdom","Pooled")



# Figure 1 (Want More or Fewer Asylum Seekers) Results --------------------

svyww$WantMoreAS <- as.numeric(svyww$AsylumHome > 0)

#stack a pooled sample
hold <- svyww
hold$cty <- "Pooled"
svywww <- rbind(svyww,hold)

#empty containers
est <- c()
se <- c()
n <- c()

for (i in 1:length(countries2)){
  
  dfrw <- subset(svywww,svywww$cty == countries2[i])
  
  surveyobject <- svydesign(ids = ~0, variables = data.frame(dfrw$WantMoreAS), data = dfrw, weights = dfrw$weight)
  est[i] <- svymean(~dfrw.WantMoreAS,surveyobject)[1]
  se[i] <- sqrt(vcov(svymean(~dfrw.WantMoreAS,surveyobject)))
  
  n[i] <- nrow(dfrw)
  
}

#Change to percentage from proportion
est <- round(est*100,2)
se <- round(se*100,2)

ddat <- data.frame(countries2,n,est,se)

alpha <- 0.05
cval <- qnorm(1-alpha/2)
ddat$lower <- round(ddat$est - cval*ddat$se,2)
ddat$upper <- round(ddat$est + cval*ddat$se,2)

#Estimate table:
ddatout <- ddat[c(order(ddat$countries2[1:15]),16),]
ddatout


# Figure 2 (Distribution of Preferences) Results --------------------------


# A. Baseline: Info Control Cons Control

svyww$PreferStatusQuo <- as.numeric(svyww$SharePref == "1. first entry")
svyww$PreferSame <- as.numeric(svyww$SharePref == "2. same for all")

svywww <- subset(svyww,svyww$InfoTrt == 0 & svyww$ConsTrt == 0)
tcondition <- "A. Baseline"

hold <- svywww
hold$cty <- "Pooled"
svywww <- rbind(svywww,hold)

est.prop <- c()
est.same <- c()
est.SQ <- c()
n <- c()

for (i in 1:length(countries2)){
  
  dfrw <- subset(svywww,svywww$cty == countries2[i])
  
  surveyobject <- svydesign(ids = ~0, variables = data.frame(dfrw$PreferProp,dfrw$PreferSame,dfrw$PreferStatusQuo), data = dfrw, weights = dfrw$weight)
  est.prop[i] <- svymean(~dfrw.PreferProp,surveyobject)[1]
  est.same[i] <- svymean(~dfrw.PreferSame,surveyobject)[1]
  est.SQ[i] <- svymean(~dfrw.PreferStatusQuo,surveyobject)[1]
  
  n[i] <- nrow(dfrw)
  
}

thisdf <- data.frame(countries2,n,est.prop,est.same,est.SQ,Condition=tcondition)
thisdfA <- thisdf


# B. Info Treatment only: Info Treated Cons Control

svyww$PreferStatusQuo <- as.numeric(svyww$SharePref == "1. first entry")
svyww$PreferSame <- as.numeric(svyww$SharePref == "2. same for all")

svywww <- subset(svyww,svyww$InfoTrt == 1 & svyww$ConsTrt == 0)
tcondition <- "B. Information Treatment"

hold <- svywww
hold$cty <- "Pooled"
svywww <- rbind(svywww,hold)

est.prop <- c()
est.same <- c()
est.SQ <- c()
n <- c()

for (i in 1:length(countries2)){
  
  dfrw <- subset(svywww,svywww$cty == countries2[i])
  
  surveyobject <- svydesign(ids = ~0, variables = data.frame(dfrw$PreferProp,dfrw$PreferSame,dfrw$PreferStatusQuo), data = dfrw, weights = dfrw$weight)
  est.prop[i] <- svymean(~dfrw.PreferProp,surveyobject)[1]
  est.same[i] <- svymean(~dfrw.PreferSame,surveyobject)[1]
  est.SQ[i] <- svymean(~dfrw.PreferStatusQuo,surveyobject)[1]
  
  n[i] <- nrow(dfrw)
  
}


thisdf <- data.frame(countries2,n,est.prop,est.same,est.SQ,Condition=tcondition)
thisdfB <- thisdf


# C. Consequences Treatment Only: Info Control Cons Treated

svyww$PreferStatusQuo <- as.numeric(svyww$SharePref == "1. first entry")
svyww$PreferSame <- as.numeric(svyww$SharePref == "2. same for all")

svywww <- subset(svyww,svyww$InfoTrt == 0 & svyww$ConsTrt == 1)
tcondition <- "C. Consequences Treatment"

hold <- svywww
hold$cty <- "Pooled"
svywww <- rbind(svywww,hold)

est.prop <- c()
est.same <- c()
est.SQ <- c()
n <- c()

for (i in 1:length(countries2)){
  
  dfrw <- subset(svywww,svywww$cty == countries2[i])
  
  surveyobject <- svydesign(ids = ~0, variables = data.frame(dfrw$PreferProp,dfrw$PreferSame,dfrw$PreferStatusQuo), data = dfrw, weights = dfrw$weight)
  est.prop[i] <- svymean(~dfrw.PreferProp,surveyobject)[1]
  est.same[i] <- svymean(~dfrw.PreferSame,surveyobject)[1]
  est.SQ[i] <- svymean(~dfrw.PreferStatusQuo,surveyobject)[1]
  
  n[i] <- nrow(dfrw)
  
}

thisdf <- data.frame(countries2,n,est.prop,est.same,est.SQ,Condition=tcondition)
thisdfC <- thisdf


# D. Both treatments: Info Treated Cons Treated

svyww$PreferStatusQuo <- as.numeric(svyww$SharePref == "1. first entry")
svyww$PreferSame <- as.numeric(svyww$SharePref == "2. same for all")

svywww <- subset(svyww,svyww$InfoTrt == 1 & svyww$ConsTrt == 1)
tcondition <- "D. Both Treatments"

hold <- svywww
hold$cty <- "Pooled"
svywww <- rbind(svywww,hold)

est.prop <- c()
est.same <- c()
est.SQ <- c()
n <- c()

for (i in 1:length(countries2)){
  
  dfrw <- subset(svywww,svywww$cty == countries2[i])
  
  surveyobject <- svydesign(ids = ~0, variables = data.frame(dfrw$PreferProp,dfrw$PreferSame,dfrw$PreferStatusQuo), data = dfrw, weights = dfrw$weight)
  est.prop[i] <- svymean(~dfrw.PreferProp,surveyobject)[1]
  est.same[i] <- svymean(~dfrw.PreferSame,surveyobject)[1]
  est.SQ[i] <- svymean(~dfrw.PreferStatusQuo,surveyobject)[1]
  
  n[i] <- nrow(dfrw)
  
}


thisdf <- data.frame(countries2,n,est.prop,est.same,est.SQ,Condition=tcondition)
thisdfD <- thisdf


#Results:
thisdfA
thisdfB
thisdfC
thisdfD



# Figure 3 (Delta in Support for Prop vs. Status Quo) Results ----------------------------------------------------------------


#First for respondents not shown consequences (pooling over info treatment)
svywww <- subset(svyww,svyww$ConsTrt == 0)

hold <- svywww
hold$cty <- "Pooled"
svywww <- rbind(svywww,hold)

est0 <- c()
se0 <- c()

for (i in 1:length(countries2)){

  dfrw <- subset(svywww,svywww$cty == countries2[i])
  
  nw <- nrow(dfrw)
  dfrw$weight <- dfrw$weight/mean(dfrw$weight)
  
  est0A <- (1/nw)*sum(dfrw$PreferProp*dfrw$weight)
  
  dfrw$WeightNorm <- dfrw$weight/sum(dfrw$weight)
  vhat.est0A <- sum(((dfrw$WeightNorm)^2) * (dfrw$PreferProp - est0A)^2)
  
  dfrw$SQPref <- as.numeric(dfrw$SharePref == "1. first entry")
  
  
  est0B <- (1/nw)*sum(dfrw$SQPref*dfrw$weight)
  vhat.est0B <- sum(((dfrw$WeightNorm)^2) * (dfrw$SQPref - est0B)^2)
  
  cov.est0AB <- sum(((dfrw$WeightNorm)^2) * (dfrw$SQPref - est0B)*(dfrw$PreferProp - est0A))
  
  est0[i] <- est0A - est0B
  
  var0 <- vhat.est0A + vhat.est0B - 2*cov.est0AB
  se0[i] <- sqrt(var0)
  
}

est0 <- round(est0*100,2)
se0 <- round(se0*100,2)


#Second for respondents shown consequences (pooling over info treatment)
svywww <- subset(svyww,svyww$ConsTrt == 1)

hold <- svywww
hold$cty <- "Pooled"
svywww <- rbind(svywww,hold)

est1 <- c()
se1 <- c()

for (i in 1:length(countries2)){
  
  dfrw <- subset(svywww,svywww$cty == countries2[i])
  
  nw <- nrow(dfrw)
  dfrw$weight <- dfrw$weight/mean(dfrw$weight)
  
  est1A <- (1/nw)*sum(dfrw$PreferProp*dfrw$weight)
  
  dfrw$WeightNorm <- dfrw$weight/sum(dfrw$weight)
  vhat.est1A <- sum(((dfrw$WeightNorm)^2) * (dfrw$PreferProp - est1A)^2)
  
  dfrw$SQPref <- as.numeric(dfrw$SharePref == "1. first entry")
  
  est1B <- (1/nw)*sum(dfrw$SQPref*dfrw$weight)
  vhat.est1B <- sum(((dfrw$WeightNorm)^2) * (dfrw$SQPref - est1B)^2)
  
  cov.est1AB <- sum(((dfrw$WeightNorm)^2) * (dfrw$SQPref - est1B)*(dfrw$PreferProp - est1A))
  
  est1[i] <- est1A - est1B
  
  var1 <- vhat.est1A + vhat.est1B - 2*cov.est1AB
  se1[i] <- sqrt(var1)
  
}

est1 <- round(est1*100,2)
se1 <- round(se1*100,2)


ddat <- data.frame(countries2=rep(countries2,2),givenshares=rep(c("No","Yes"),each=16),est=c(est0,est1),se=c(se0,se1))

ddat$lower <- ddat$est - 1.96*ddat$se
ddat$upper <- ddat$est + 1.96*ddat$se

#Results

ddat

#pvalues for differences:

2*(1-pnorm(abs((est1-est0)/sqrt(se0^2 + se1^2))))

